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Ultrarelativistic collisions of black holes are ideal gedanken experiments to study the nonlinearities 
of general relativity. In this paper we use semianalytical tools to better understand the nature of 
these collisions and the emitted gravitational radiation. We explain many features of the energy 
spectra extracted from numerical relativity simulations using two complementary semianalytical 
calculations. In the first calculation we estimate the radiation by a "zero-frequency limit" analysis 
of the collision of two point particles with finite impact parameter. In the second calculation 
we replace one of the black holes by a point particle plunging with arbitrary energy and impact 
parameter into a Schwarzschild black hole, and we explore the multipolar structure of the radiation 
paying particular attention to the near-critical regime. We also use a geodesic analogy to provide 
qualitative estimates of the dependence of the scattering threshold on the black hole spin and on 
the dimensionality of the spacetime. 
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I. INTRODUCTION 

Solving Einstein's equations numerically is a highly 
nontrivial task. After the recent breakthroughs in numer- 
ical relativity (NR) binary black hole (BH) mergers 
have been routinely carried out by several groups world- 
wide. As is often the case in physics, the results of nu- 
merical simulations are best understood or interpreted 
by studying simplified models that capture the main fea- 
tures of the problem. In this paper we will argue that 
semianalytical tools are particularly useful to gain in- 
sight into the fascinating but complex problem of ultra- 
relativistic BH collisions. 

High-energy BH encounters are a difficult undertaking 
m NR (see Refs. 0-0] for a discussion of ultrarelativis- 
tic collisions in four dimensions, Refs. fcllllT] for recent 
progress on D-dimensional simulations and Ref. (l2| for 
a review of the most outstanding questions in ultrarela- 
tivistic BH collisions before the recent breakthroughs in 
NR) . Efficient adaptive mesh refinement and wave extrac- 
tion techniques are required because the problem involves 
various scales: the BHs are "pancake-shaped" because 
of Lorentz contraction, the large speeds involved in the 
collision require large initial separations to define asymp- 
totic states, and high resolution is required to study the 
dynamics of the final BH (if any) formed as a result of 
the merger. Further difficulties arise from the spurious 
radiation present in the initial data. 

Despite their relatively limited accuracy, simulations 
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of ultrarelativistic BH collisions to date have provided 
definite predictions for nonspinning BHs in four dimen- 
sions 0t3- For example these simulations show that, 
even in the highly symmetric case of head-on collisions, 
as much as ~ 14 ± 3% of the energy of the system can 
be radiated in gravitational waves [H, Q. They further 
demonstrate the existence of three distinct regimes de- 
pending on the impact parameter b: immediate mergers, 
nonprompt mergers and the scattering regime. These 
regimes are separated by two special values of the im- 
pact parameter: the threshold of immediate merger b* 
and the scattering threshold 6 sca t. Roughly speaking, for 
b < b* merger occurs within the first encounter, whereas 
for b* < b < b scat it docs not, but sufficient energy is 
radiated to put the binary into a bound state that even- 
tually results in a merger. For the largest initial center 
of mass velocities of v = 0.94 studied to date, as much as 
~ 35 ± 5% of the energy can be carried away by gravita- 
tional radiation 0, Q . 

Close to the threshold of immediate merger, the bi- 
nary exhibits "zoom-whirl" behavior, an extreme version 
of relativistic perihelion precession which can be pre cisely 
defined and understood in the geodesic limit [IJ-QjI . For 
point particles orbiting BHs, zoom-whirl orbits are inti- 
mately related to the existence of unstable spherical or- 
bits (unstable circular orbits at radii 3M < r < 6M 
in the special case of Schwarzschild BHs). The very ex- 
istence of zoom-whirl orbits in comparable-mass BH en- 
counters is perhaps a strong hint that simple, semianalyt- 
ical approaches (in this case, the study of point-particle 
geodesies around BHs) can provide valuable insight into 
the general solution of the problem. 

NR simulations of ultrarelativistic BH collisions have 
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FIG. 1: Energy spectrum for the dominant (quadrupolar, i.e. 
I — 2) component of the gravitational radiation computed 
from NR simulations of the head-on collision of two equal- 
mass BHs (from Q). The collision speed in the center-of- 
mass frame, /3 = v/c, is indicated in the legend. The en- 
ergy spectrum is roughly flat (independent of frequency) up 
to the quasinormal mode (QNM) frequencies (marked by ver- 
tical lines), after which it decays exponentially. All quantities 
are normalized to the Arnowitt-Deser-Misner (ADM) mass of 
the system Madm- 



also provided information on the structure of the Fourier- 
domain energy spectra dE lad / 'duj . An interesting feature 
of the spectra of comparable-mass BH collisions is the 
presence of a nonvanishing zero-frequency limit (ZFL). 
For head-on collisions, Fig. [T] shows that the energy spec- 
trum is roughly flat (independent of frequency) up to 
some cutoff frequency, after which it decays exponen- 
tially (in this figure and elsewhere in the paper we use 
geometrical units G = c = 1). Grazing ultrarclativistic 
BH collisions also lead to energy spectra with complex 
features, although this was not discussed in Ref. @ be- 
cause of space limitations. These complex spectra are 
discussed further in section [V] of the present work. 

One of the motivations of this paper is to explore sim- 
plified, analytical and semianalytical descriptions of ul- 
trarelativistic BH mergers that provide qualitative (and 
sometimes even quantitative) explanations of the main 
features of these spectra [||. In particular, we shall con- 
centrate on the combination of the following two ap- 
proaches: a linearized calculation in flat spacetime pi- 
oneered by Weinberg and Smarr [l8l - [2l| and BH pertur- 
bation theory |23 - l28j . We will refer to the Weinberg- 
Smarr approach (which is valid for arbitrary velocities, 
as long as the radiated energies are small) as the "ZFL 
calculation," because it provides a good approximation 
of the emitted energy spectrum at low frequencies. The 
ZFL calculation assumes an instantaneous collision in flat 



spacetime, and in four dimensions it yields a flat spec- 
trum at all frequencies - precisely what one would ex pec t 
from the analogous problem in electromagnetism [29j . 
The origin of the exponential cutoff in the NR spectra 
can be understood through the second approach, i.e. by 
computing the radiation from point particles falling in 
a curved BH spacetime [22H281 ] . These perturbative cal- 
culations show that the radiation in a given multipolc 
"shuts off" at the real part of the lowest QNM frequency 
of the BH for the multipole in question [30, UK ■ 

In this paper we address several questions related to 
the puzzling behavior of ultrarclativistic BH collisions: 

1. We use geodesic calculations to explore the de- 
pendence of the scattering threshold on BH spin 
and dimensionality of the spacetime. This thresh- 
old determines whether the two BHs will eventu- 
ally merge or not, and it is of fundamental im- 
portance to estimate cross sections for BH produc- 
tion in high-energy collisions [32j . Classical general 
relativity in D dimensions should be adequate to 
determine the cross section for BH production at 
trans-Planckian collision energies and the fractions 
of the collision energy and angular momentum lost 
in gravitational radiation. This information will be 
of paramount importance to improve the model- 
ing of microscopic BH production in event genera- 
tors such as Truenoir, Charybdis2, Catfish or 
BLACKMAX [3314361 ] . The event generators will then 
provide a description of the corresponding evapo- 
ration phase, which might be observed during LHC 
collisions. 

2. We reproduce the qualitative (and sometimes quan- 
titative) features of the energy spectra by a "bare- 
bones" approach, where we replace the colliding 
BHs by two point particles and we estimate the 
radiation by a ZFL calculation with finite impact 
parameter. The generalization to finite impact 
parameter requires a specific model to bind the 
particles after the collision and guarantee energy- 
momentum conservation 1371 ] . but our qualitative 



conclusions should be roughly independent of the 
details of the model. 

3. We gain insight into the details of the emitted radi- 
ation (particularly in the near-critical regime) by a 
perturbative study of the ultrarclativistic infall of 
point particles with generic impact parameter into 
a Schwarzschild BH. 

The two-pronged analytical approach above provides 
a better understanding of most of the main features of 
the energy spectra. In particular, we reach the following 
conclusions: 

i. In both the extreme-mass ratio case and in 
comparable-mass NR simulations, the ZFL is by and 
large independent of the impact parameter b. 
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ii. For a given multipolar index I, the slope of the differ- 
ent m-components at low frequencies is proportional 
to b. 

iii. The exponential cutoff is related to the QNM fre- 
quencies of the final BH. 



threshold in the Schwarzschild background. Then we pro- 
ceed to the computation of the critical impact parameter 
in a four-dimensional rotating (Kerr) background, and 
finally we generalize the analysis to higher-dimensional 
rotating (Myers-Perry) backgrounds. 



iv. Higher multipoles of the radiation become increas- 
ingly important as the binary becomes ultrarelativis- 
tic, leading from an exponential scaling of the form 
e ~a R i j n .j. ne nonre i a tivistic limit to a scaling of the 
form bn/l in the ultrarelativistic limit, where (an, bn) 
are constants. 

v. The total radiation is proportional to the number of 
orbits near the scattering threshold for nonreiativistic 
infalls, but it is enormously enhanced (essentially by 
a resonance with the QNM frequencies of the final 
BH) when the collision is ultrarelativistic. 

The plan of the paper is as follows. In section [H] 
we compute the critical impact parameter for geodesies 
in four and higher dimensions. In section Mil we ex- 
plore a simple toy model to estimate the radiation in 
the ZFL for collisions with finite impact parameter. In 
section IIVI we compute the energy, angular momentum 
and linear momentum radiated by point particles falling 
into Schwarzschild BHs in four dimensions with arbitrary 
energy. In section [V] we show how the ZFL and point- 
particle calculations shed light on some features of NR 
simulations of ultrarelativistic BH collisions. We con- 
clude by pointing out possible extensions of our inves- 
tigation. Appendix [X] gives some technical details on 
the Sasaki-Nakamura formalism, and Appendix [B] shows 
how the ZFL results can be decomposed in terms of spin- 
weighted spherical harmonics. 

We work mostly in four spacetime dimensions with 
signature (—,+,+,+) [Hj], except for the discussion in 
section III C[ where we explore higher-dimensional BHs. 
Greek letters (fx, v, . . .) in index lists range over all space- 
time indices. The Einstein summation convention is em- 
ployed unless otherwise specified. Overhead dots stand 
for partial differentiation with respect to the proper time 
t, E = dE/dr; 0(A) stands for terms of order A. 



II. SCATTERING THRESHOLD IN BLACK 
HOLE SPACETIMES 

One of the main ingredients to estimate BH produc- 
tion rates in particle accelerators is the cross section for 
BH production (32j . Several analytical ap prox imations 
have been used to estimate cross sections [T2| - In this 
section we explore a particularly simple approach. We 
compute the critical impact parameter between plunge 
and scattering for particles following geodesies in dif- 
ferent BH metrics, in order to estimate the qualitative 
dependence of the cross section on BH spin and dimen- 
sionality. We begin with a brief review of the scattering 



A. Schwarzschild black holes 

Geodesies in a Schwarzschild background are com- 
pletely determined by their energy and (z-component of) 
orbital angular momentum per unit rest mass. It is useful 
to replace the total energy E and the orbital angular mo- 
mentum L z by the parameters E = Ejfx and L z = L z / fx. 
Here ix is the particle's rest mass and E is related to the 
velocity v of the particle at infinity by E = (1 — v 2 )^ 1 ^ 2 . 
The geodesic radial behavior is governed by the rela- 
tion (dr/dr) 2 = E 2 — V e g(r, L z ), where r is the proper 
time and the effective potential for a Schwarzschild back- 
ground is Vcs = (1 + L 2 Jr 2 )f, with / = 1 - 2M/r and 
M the BH mass. 

Geodesies can be classified according to how their en- 
ergy compares to the maximum value of the effective po- 
tential, which in this case is given by 



ymax 
en 




The scattering threshold is then defined by the condition 
E 2 = y c 'g ax : unbound orbits with E 2 > V^g* are cap- 
tured, while those with E 2 < V^g ax are scattered. Given 
some E, the critical radius or impact parameter & cr ; t that 
defines the scattering threshold is obtained by solving the 
condition E 2 = V^ &x for L z = L cl - it , and then using the 
relation 



frcrit — L cr it(E 2 — 1) 



-1/2 



(2) 



Note that we mostly work in the point particle approxi- 
mation with no radiation reaction effects, so we adopt a 
slightly different terminology from that of Ref. @ : here 
L C rit denotes the critical angular momentum separating 
plunging trajectories from scattering trajectories. 

Ultrarelativistic collisions can be modeled by large- 
energy geodesies. One can then show that a high-energy 
orbit plunges when 



^crit 
M 2 M 2 



^-9-^ + 0(1/^), (3) 



where the last relation is valid in the limit E 2 3> 1. The 
condition for scattering is L z > L cr it . One can also easily 
show that all orbits plunge if L 2 < 12M 2 when the par- 
ticle is at rest at infinity (E = 1). If E > 1, scattering 
orbits exist only if L z > AM . 
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B. Kerr black holes 

The radial motion of equatorial geodesies in a Kerr 
background can also be described in terms of an effective 
potential [3^. In this case, however, the potential for 
timelike geodesies takes a more complicated form, which 
we can parametrize as 



j ( T 1 r 3 ) ( ) 



The Kerr spin parameter a and the reduced spin param- 
eter j are related to the spin angular momentum J via 
J = Ma = jM 2 . As before, E and L z stand for the en- 
ergy and (z-component of) angular momentum per unit 
rest mass. For convenience, we have also defined the 
energy-dependent "spin-deformation" parameters 



o-i = AM 2 E , 



M (1 - E 



03 



-2M 3 E 2 



(5) 

For fixed j, orbits with a turning point are identified 
by pairs (E, L z ) for which V e g has a double root. The 
extrema of the effective potential are defined by 



dV eS _ 2M + 6ML 2 Z 2L 2 Z 3ja 1 L z . 2 /2<j 2 + 3cr 3 



and correspond to radii 
L 2 



(6) 



2M 



J L 2 2M 



. 2 ,'2(72 . 6M<r 3 
- 3 I JT + —FI- 

±J -, J-J ~ 



J Li 



12M 7 

1/2 



+ 3 



&Ma x 



(7) 



L cr it, and then using Eq. ([2j . Figure [2] plots this critical 
impact parameter as a function of energy for different 
values of the spin, where we have solved for 6 cr it numeri- 
cally. Observe that & cr it asymptotes rapidly to a constant 
as the energy approaches E ~ 2. For large, positive 
spins (co- rotating case), the critical impact parameter 
asymptotes to 6 cr it ~ 2.33M for E = 2, while for large, 
negative spins (counter-rotating case) it asymptotes to 
&crit ~ 7.62M. As the spin is decreased, the critical im- 
pact parameter for corotating geodesies increases, so that 
the smallest impact parameter corresponds to maximally 
spinning BHs. 

The above considerations suggest that an asymptotic 
analysis might allow for a better analytical understanding 
of the critical impact parameter in the ultrarelativistic 
regime. Let us then define the perturbation parameters 

1 > e = 1 - j > , 1 > £ = 1/E > . (8) 



Performing an asymptotic expansion to second order in 
e, the solution to E 2 = V^J ax becomes 

2 , 2 (C 2 - 1) 



icrit ~ ^ + V e_ 



e (c 3 - 3) 



(9) 



= 



e 3/2 



(e-i) (3e 4 + 5)v^2F +0(e2) 



€(e 2 -3) j 

Rccxpanding this equation to second order in £, we find 

2 5^ 
— e H 

3 108 



(10) 



1 + ^_J e+ ^. W 0(e >,a 



The maximum of the effective potential is then equal to 



T/max 
'eff 



14ff(r = r_), which generalizes Eq. ([T]). 




Eq. (|10[) is formally a bivariate expansion, i.e. an ex- 
pansion in two independent perturbation parameters: e 
and £. This expansion is a fractional Frobenius series, 
with the regular limit 6 cr i t — > 2 when e — > and ^ — > 0. 
The energy parameter contributes only to second order, 
while the spin contributes at fractional leading order. A 
comparison of this estimate to the numerical solution of 
Fig.[2]shows good agreement: when j = 0.999 and E = 2 
Eq. (|T0|) predicts 6 cr ; t ~ 2.335M, while the numerical re- 
sult is &num = 2.396Af. Therefore the relative fractional 
error of the asymptotic expansion (dominated by the ne- 
glected relative E~ 3 terms) is approximately 2.6%. 

In the ultrarelativistic limit (£ — > 0), we can solve the 
critical impact parameter equation exactly for arbitrary 
j to find [39j 



frcrit . , _ 

— — = —J + o COS 
M 



arccos(— j) 



(11) 



FIG. 2: Critical impact parameter versus reduced energy E 
for different values of j in a Kerr background. 

If we are given some values of (E, j), the critical impact 
parameter & cr i t can be obtained by solving E 2 = V^ ax for 



or alternatively b cr \t/M = —j + 3x + 3/x, where x = 
[—j + (j 2 — l) 1 / 2 ] 1 / 3 . The asymptotic expansion is very 
accurate in the ultrarelativistic limit: for example, for 
j = 0.999 this exact formula gives b clit /M = 2.07812987, 
to be compared with b clit /M - 2.07812992 from Eq. (JTD]). 
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C. Myers-Perry black holes 

We now generalize the above analysis to higher- 
dimensional BHs, considering for illustration the Myers- 
Perry solution with a single angular momentum direction 
[4pj . The metric of a Z?-dimensional Kerr BH with only 
one nonzero angular momentum parameter is given in 
Boyer-Lindquist-type coordinates by [40[ 



, , A - a 2 sin 2 , 9 2a(r 2 + a 2 - A) sin 2 d 
ds A = dr 



(r 2 + a 2 ) 2 - Aa 2 sin 2 i? 



sin 2 ddtp 2 



-dr 2 



Y>dfi 2 +r 2 cos 2 $dn 2 D _ 4 , 



where 



S 
A 



r 2 + a 2 cos 2 
r 2 + a 2 - M»r 



5-D 



dtdtp 



(12) 



(13) 
(14) 



dft 2 j_ 4 denotes the standard metric of the unit (D — 4)- 
sphere, and M* is related to the BH mass. In fact this 
metric describes a rotating BH in an asymptotically flat 
vacuum spacetime with mass AA and angular momentum 
J given by 



M 



D-2 A 1 

-A {D _ 2) M*, J = — A {D _ 2) M*a, (15) 



16tt 



where Am_2) is the area of a unit {D — 2)-spherc: 



-4 



(D-2) 



r[(D-i)/2] 



(16) 



Timelikc equatorial geodesies depend on the effective 
potential 



(a 2 E 2 - L 2 ) - JiAi 



(17) 

where D is the dimensionality of spacetime and Ac = 
r 2 + a 2 - M*r 5 ~ D . For D = 4, the quantity M* and 
the mass of the BH are related via M* = 2Ai. In higher 
dimensions M* does not have units of length, but instead 
[M„] = (length) 13-3 . We also introduce a dimensionlcss 
spin parameter 



J 



(M*/2) 1 /(^-3) ' 



(18) 



which reduces to the corresponding Kerr quantity for 
D = 4. 

When D = 5 we can carry out an asymptotic analysis 
similar to the one presented in the previous subsection 
(see also [!l|). In the ultrarelativistic limit (E — >• oo) we 
have 



dD=5) 



(19) 



which is to be contrasted with the four-dimensional result 
presented in Eq. pT|) . This asymptotic analysis reveals 
that both in four and five dimensions the overall effect of 
spin for corotating geodesies is to reduce the critical im- 
pact parameter, suggesting that the highest energy emis- 
sion occurs for maximally rotating BHs. 
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FIG. 3: Critical impact parameter for E — 10 as a function of 
j for Myers-Perry BHs in different dimensions D. The inset 
zooms along the j/-axis to show the local minima more clearly. 



For dimensions higher than five, an asymptotic analy- 
sis is more involved due to the higher inverse polynomial 
order of the effective potential. However, a numerical 
calculation is straightforward. In Fig. [3] we plot the di- 
mensionless impact parameter 6 cr it/(Af*/2) 1 ^- D ~ 3 ' as a 
function of j for different values of D. The calculations 
in Fig. [3] refer to E — 10, but this is already a good 
approximation of the asymptotic value of the impact pa- 
rameter in the ultrarelativistic region. Observe that as 
D increases, the minimum in & cr it/(M*/2) 1 /(- D ~ 3 - ) corre- 
sponds to smaller values of j. 

The (v,j) phase space is explored in more detail in 
Fig. IU where we show contour plots of the dimension- 
less impact parameter for both corotating and counter- 
rotating geodesies in different dimensions. Note that the 
smallest critical impact parameter corresponds to coro- 
tating, ultrarelativistic geodesies, but the BH must not 
necessarily be maximally spinning when D > 5. For 
6 < D < 10 the critical impact parameter presents a 
local minimum 6™" as a function of j, which is well de- 
scribed by the following quadratic fit: 



3.978 - 0.679L> + 0.033D 2 . (20) 



This fit was constructed considering only spacetimes with 
D < 10, so its extrapolation to D > 10 should be used 
with caution. 
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Critical impact parameter, D=4 



Critical impact parameter, D=5 




FIG. 4: Contour plots of the critical impact parameter in the (v,j) plane for corotating orbits (j > 0) and counterrotating 
orbits (j < 0) in dimensions D = 4, D = 5, D — 6 and D = 7. 



III. THE ZFL IN GRAZING COLLISIONS 

In this section we generalize the classic ZFL calcula- 
tions for head-on collisions [13, H2 to the case of collisions 
with finite impact parameter. The initial configuration 
consists of two point particles with mass freely mov- 
ing toward each other with constant, positive velocity 
ffe, corresponding to boost factors Ek = (1 — Wfe) -1 ^ 2 
(k = 1,2). For convenience the axes are oriented such 
that the initial motion is in the x-direction (see Fig. [5J . 
We assume that at x — the particles "collide" with 
generic impact parameter b and form a single final body 
(strictly speaking this assumption is only valid for small 
impact parameters, because we expect the bodies to scat- 
ter when b is large enough). Since the collision is not 
head-on (and since the energy loss is not included in 



the motion of point particles), some confining force is 
necessary to bind the particles. In fact, we show below 
that additional "stresses" are required to guarantee en- 
ergy conservation (cf. Ref. (37|). 

Before the collision the particles have four-positions 
and four-momenta given by: 



Pi 



(Mit.fc.O), 
£iMi(Mi,0,0) 



x 2 = (*) -vit, -6 

P2 



0), (21) 



" - E 2 M 2 (l,-v 2 ,0,0). 



where £i (—£2) is the projection of the position of particle 
1 (2) along the y-axis before the collision. If the system's 
center of mass is at y = 0, the position of particle 2 before 
the collision can also be written as y = — (b — £1) = — £ 2 , 
where b is the impact parameter. 

At t = the particles become constrained to move as 
if they were attached to an infinitesimally thin, massless 



7 



t<0 



t>o 



FIG. 5: The system before and after the collision, £1 is defined 
in Eq. p2]l. 



rod of length b. This fictional rod is an idealization, but 
it is necessary to guarantee energy-momentum conserva- 
tion. For t > the particles remain attached to the rod, 
so that (in the center-of-momentum frame) they rotate 
around the origin at fixed separation b. Using primes to 
denote final states, the four-positions and four-momenta 
after the collision arc: 

= (t,^s,^c,o), 

(a£)' = (t,-&S,-bC,0), 

(p?)' = E 1 M 1 (1,60(7,-^05,0), 

(p%y = e 2 m 2 (i,-&qc, 605,o), 

where S = sin (fit), C = cos (fit). 

A few relations can be derived between (6j6) an d 
(vi,v 2 ) due to the constraints we imposed. The require- 
ment that the system's center of mass is at y = implies 
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bE 2 M 2 



(22) 



E X M X +E 2 M 2 

In the center-of-momentum frame we also have 

EiMxVi = E 2 M 2 v 2 . (23) 

For an instantaneous collision the energy-momentum 
tensor of the system is given by: 



k=l 



|^<5 3 (x-x fc (t))6(-t) (24) 



+ ^W|MM 4 . (x _ 4(i))e(() , 



fe=i 



where the boldface denotes a three- vector, and the angu- 
lar momentum is 



S 3 = J (x^ 20 - x 2 T 10 )d 3 x 



(25) 



• J2 [®(-t)E k M k v k £ k + e(t)E k M k fli 



fc=i 



Angular momentum conservation implies that the rota- 
tion frequency must be 







EiM±vi^i + E 2 M 2 v 2 £ 2 
ExM^l + E 2 M 2 J\ 



(26) 



If we naively take the stress-energy tensor of Eq. ((24)) to 
be the full energy-momentum of the system, we would 
find that it is not covariantly conserved, i.e. V = 
for v = t, z but V^T^" 7^ for v = x, y. In fact, one 
finds (in the center-of-momentum frame) 



-E 1 M 1 ^ 1 n 2 SS(x - £iS)5(y - ^C)S(z)Q(t) 
+ E 2 M 2 ^ 2 fl 2 SS(x + &S)5(y + &C)5(z)Q(t) , 
V M T^ = -E 1 Mi^fl 2 C5(x - £iS)6(y - ^C)6(z)Q(t) 
+ E 2 M 2 t; 2 n 2 C6(x + &S)5(y + t; 2 C)5(z)0(t) . 

Physically, this nonconservation of stress-energy is due 
to neglecting the energy-momentum associated with the 
fictitious rod that keeps the particles in circular orbit. 

Energy-momentum conservation can be enforced by 
adding an additional term for each particle that repre- 
sents this constraining force. The contribution of such 
forces to the gravitational radiation emitted by a parti- 
cle in circular orbit was studied by Price and Sandberg 
[37| . By adding a radial tension T k (r) for each particle 
and imposing that V M T A " y = we get the following con- 
tributions to the energy-momentum tensor: 



3S£,(t,x) = -S 2 8(co S 6)Q(t)J2T k (r)5(4> + nt-(t> k ), 

2 

7T ns (*,x) = -C 2 5(co S 6)Q(t)J2n(r)6(<f> + nt-<f> k ), 

fc=i 
2 

TZ s (t,x) = -SCS(cos9)e(t)Y,r k (r)6(^ + flt-^ k ), 



k=l 



where </>i = tt/2, 4> 2 — 3n/2 and 
M fc ye(ft-r) 



Here r = ^ x 2 + y 2 + z 2 , 8 is the polar angle measured 
from the positive z-axis, and <f> is the azimuthal angle 
in the x-y plane measured from the £-axis (see Fig. 
The factor [l — (£fcO) 2 ] is just the boost factor for 
particle k in circular motion with angular frequency O. 
The stresses vanish for b = 0, as one would expect. 



A. Equal-mass collisions 

In this subsection we study the equal-mass case M /2 = 
Mi = M 2 , where M is the total mass, and V\ = v 2 = v, 
E\ = E 2 = E. According to Eq. ([26]) . after the collision 
the particles are on a bound circular orbit with radius 
6/2 and rotational frequency O = 2v/b. The Fourier 
transform of the energy-momentum tensor (|24[) yields 
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T^Ck UJ) = P ^ P " p -»fcyf»/2 j P2P2 ik v b/2 , / 28 N 

V ' ' 2-kiExM^uj - v x k x ) 2ttiE 2 M 2 (uj + v 2 k x ) K ' 



where d 4 x = dx dy dz dt and k is the wave vector: 



Radiation Spectrum 



k x = lu sin <p cos 9 , k y = uj sin </> sin 9 , k z = uj cos < 



We also have 



w6 



e exp iAfc — sin 9 sin(Qi 



(29) 



(30) 



where fc = 1, 2 is the particle index, Ai = — 1 and A2 = 
1. If we set a = fit + <f> and 277^ = X^Lub sin 9 the last 
exponential can be written in terms of Bessel functions 
of the first kind, using the Jacobi- Anger expansion (42j : 



n— +00 



(31) 



n— — 00 



For large n the Bessel functions satisfy \A1 



Jn{v) 



(32) 



A time- integration introduces an additional factor of 1 /n, 
so the series converges rapidly for large \n\ and we can 
truncate it at some moderately large value of n = N 
to get an accurate approximation of the integral 1 . Typi- 
cally, N > 10 is sufficient for an accuracy of 1% or better. 

The integration of the stresses proceeds in a similar 
way. After integrating in 9 and </>, the same Bessel func- 
tion expansion can be used for the time-integration. The 
integral of Bessel functions with respect to r can be eval- 
uated using the following identity [42| : 



Jv{r)dr 



2 E J - 

fc=0 



2fe+i(r?), Re(i/)>-l. (33) 



1 Actually, the series should be approximated by summing from 
n = no — N to n = no + N, where no is the value of n which 
maximizes the absolute value of the terms being summed. After 
the integration terms of the form X/(u!—nO.) appear. This means 
that the largest contribution to the sum corresponds to some 
no 7^ 0. However it can be checked that N 2> no for the range 
of parameters considered here, so the sum can be taken in a 
symmetric interval around 0. 



The energy per solid angle and per unit frequency emit- 
ted in the direction k = k/w is [43| : 



d 2 E 
dudfl 



1 



= 2u/ T""(k,a,)T* (k,w) - - T\(k,w) 



(34) 

where the asterisk stands for complex conjugation. The 
energy can also be expressed in terms of the purely space- 
like components of T^ v . The conservation equation for 
T» v implies that k^T^{k,uj) = 0, so it is possible to 
write T 00 and T° l in terms of T ij : 



T 00 (k,iu) = fc i ^T ii (k,w) I 
T w (k,a>) = -fr'Tyfow). 



With these identities at hand, Eq. 
as: 



(35) 
(36) 



can be written 



d 2 E 
diodfl 



2uj 2 K ljlm (fc)T*« (k, uj)T lm (k, «) , (37) 



where following Ref. [43j we defined 



^ijlm(k) = S t lS 



1 QikjkniS^i H~ ^kihjkik-n 



ij "Im 



(38) 



In Fig. [5] we plot the energy spectrum for E = 3 along 
four different directions: k = e x , e y , e z , (e x + e y + 
v / 2e z )/2. For dimensional reasons, there is no need to 
fix b as long as the energy is plotted function of 

uj/fl. All spectra diverge when w = 2J7, as expected of 
a rigid symmetric body rotating with angular frequency 
ft. For k = e z the spectrum only diverges at oj = 2f2 
(see Ref. [44[ for a discussion of particles in circular orbit 
in the Schwarzschild geometry), but in all other direc- 
tions the spectrum diverges at even multiples of the ro- 
tational frequency Q. The same qualitative features hold 
for higher boost parameters. 



2. Head-on collisions 

Let us consider the b = limit, which corresponds to a 
head-on collision and for which we can compare against 
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(e x +e +2 e z )/2 



FIG. 6: Energy per unit solid angle and per unit frequency 
emitted in the directions k = e x , e y , e z , (e x + e y + v2e z )/2 
by equal-mass binaries with E = 3, as a function of uj/Q,. 



known results [20ll2l|. In this limit, the only nonvanish- 
ing components of the energy-momentum tensor are: 



27rwT"(k,w) = iEM 
2irujT tx (k,uj) = - 
2irwT xx (k,u) = - 



iEM 



1 — v 2 sin 2 9 cos 2 </> 
iEMv 2 sin#cos</> 
1 — v 2 sin 2 9 cos 2 cj> ' 
iEMv 2 



1 — v 2 sin 2 9 cos 2 



The energy spectrum per unit solid angle is then given 
by 



d 2 E E 2 M 2 v A (sin 2 6»cos 2 0- l)' 



dujdfl 47r 2 (■ v 2 sin 2 e cog 2 



(39) 



This agrees with previous results in the literature. In- 
deed, after a trivial redefinition of angles, Eq. (I39|> is 
equal to Eqs. (2.19) and (2T2) of Refs. [20[ andgl, 
respectively. In Refs. 0, [2l| the angular variable 9 is 
the angle between the radiation direction and the mo- 
menta of the particles; the substitution of their cos 9 by 
sin6*cos(/) yields Eq. (|3"9")l . For ease of comparison with 
NR results [1, Q , it is convenient to expand the above ex- 
pression in spin-weighted spherical harmonics with spin 
weight s = —2. Such an expansion is discussed in Ap- 
pendix IbI 

Recalling that e x corresponds to (6 = w/2 , (f> = 0), e y 
corresponds to (9 = n/2,(f) = tt/2) and e z corresponds 
to 9 = we get: 



d 2 E 
dujdQ 

d 2 E 
dujdQ 



along e x , 

E 2 M 2 v 4 



Air 2 



along e y , e z 



(40a) 
(40b) 



The radiated momentum per unit frequency for this 
head-on collision is found by the following integral over 



a two- sphere at infinity, S a 
origin: 

m _ 

dio 



centered on the coordinate 



d 2 E 
dwdQ 



riidtt , 



(41) 



where is a unit radial vector normal to Soo. In the 
present case we find that the radiated momentum van- 
ishes, i.e. dPijduj = 0, as one would expect. 



3. Zero-frequency limit 

For arbitrary impact parameters our results show that, 
in the limit boj — > 0, the energy spectrum is independent 
of b and given by Eq. (f3"5]) . This is of course consistent 
with the head-on results of Smarr [2(| and Adler and Zeks 
[2l| . Numerical calculations support this conclusion and 
reveal additional details for small but nonzero frequen- 
cies. The stress terms give the following contributions to 
the energy-momentum tensor: 

ib 2 EMn 2 



w=0 



^tens(k,w) 
wT tcns( k > W ) 



ui=0 



16tt 



= 0. 



(42) 



For lu = the constraining forces provide a nonvanish- 
ing contribution to the energy-momentum tensor. It is 
this particular contribution that allows one to recover the 
ZFL of the energy spectrum, Eq. ([39]) , for any impact pa- 
rameter. This is one of the most intriguing results of this 
incursion into the properties of the ZFL for collisions with 
nonzero impact parameter. 
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Mm 
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FIG. 7: Normalized energy spectrum per unit solid angle 
emitted in the directions k = e y (thick lines) and k = e z 
(thin lines) as a function of Mlu for several values of b/M (as 
indicated in the legend) and E — 3. As indicated by Eq. (|40p 
the ZFL for these two different directions is the same and 
approximately equal to 0.18013M 2 . 

The spectra for small frequencies along the directions 
k = e y and k = e z are plotted in Fig. [7] for different val- 
ues of b/M. All spectra have the same ZFL, as discussed 



10 



above. Specializing Eq. (|40b[) to the case E = 3 we get 
S^L = 0.18013M 2 for e y ,e z . This is in very good 

uj— 

agreement with the numerical results shown in Fig. [7J 
For small but finite frequencies, we find that the slope 



of the energy spectrum depends on direction. It is pos- 
itive (negative) for k = e y (e z , respectively), and it in- 
creases with b/M. An expansion of the energy for small 
bto yields 



d 2 E 
diodfl 



E 2 M 2 v 2 



16-k 2 (l — v 2 sin 2 9 cos 2 
1 r 



4u 2 (sm 2 (9cos 2 0- 1) 



v 2 (Sv 2 + 3) sin 6 9 cos 4 + sin 2 6 ((8v 2 + 3) cos 20 + 12 (v 2 + l)) - 
in 4 9 cos 2 (6w 4 + (2v 2 + 3) v 2 cos 20 + 20v 2 + 3) - 12 (bw) 2 j +0 [(6w) 3 ] 

I 



6 

sin 



(43) 



Thus, within our model the spectrum typically has 
quadratic corrections, except along k = e x , in which case 
the first nonvanishing contribution to the energy spec- 
trum is of order (fow) 4 . In fact, we find 



d 2 E 



k=e x 



E 2 M 2 (3 - 2v 2 ) 
9216tt 2 



■{bujf + O [{buf 



dujdQ, 

(44) 

These analytical expressions are in good agreement with 
the results shown in Fig. [7] 



B. Extreme- mass ratio collisions 

We now study collisions for fi = Mi <C M2 = M 
(the qualitative features of the radiation for generic mass 
ratio are very similar to the extreme- mass ratio case). 
The energy spectrum can be computed in the center-of- 
momentum frame. Since particle 2 is much heavier than 
particle 1, particle 2 is practically at rest in this frame, 
although we shall not neglect its motion when we com- 
pute the energy-momentum tensor. Therefore, we let 
v\ = v 3> v-2 and E\ = E. 

From Eqs. (f2"2"j) and (|26p . the angular frequency and 
position of particle 1 are given by 



E 1 ^ 1 +E 2 M^ 2 



6 



bE 2 M 



E U i + E 2 M 



(45) 



Once again we must add the stresses needed to constrain 
the particles in their orbits, in order to have a conserved 
energy-momentum tensor. 



1. Radiation Spectrum 

We compute the radiated energy using Eq. (j34|) or 
Eq. (|37|). We expand the energy-momentum tensor in 
powers of n/M, and compute the energy keeping only 
leading-order contributions in fi/M. A calculation of the 



radiation for E = 3 (in the center-of-momentum frame) 
along several different directions yields the spectra shown 
in Fig. |U 




0.2 0.4 0.6 



1 1.2 1.4 1.6 1. 

(£>IQ 



FIG. 8: Normalized energy spectrum per solid angle emit- 
ted in the directions k = e x , e y , e z , (e x + e y + v^2e z )/2 as a 
function of uj/Q. in the extreme- mass ratio case. 

The extreme-mass ratio configuration loses the angular 
symmetry of the equal-mass case. Therefore, the spectra 
now diverge for all multiples of fl for k = e x , k = e y . The 
behavior is similar for other directions. For k — e z , the 
spectrum only diverges for uj = 2Q, in agreement with 
Poisson's findings for particles in circular orbit around 
BHs 0. 



2. Head-on collisions 

For extreme-mass ratio head-on collisions {b = 0) we 
find 



d 2 E E 2 fi 2 v 4 (sin 2 cos 2 0-1) 



dujdQ 



4:Tt 2 (v sin 6 cos — l) 2 



(46) 
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This expression coincides, as it should, with Eq. (2.17) 
of Smarr [2(j, once we take into account the different 
convention on angles by an appropriate redefinition of 
angular variables (cf. the discussion following Eq. (1391) . ) 
For ease of comparison with perturbative results of point 
particles in BH spacetimes, in Appendix [Bl we compute 
analytically the multipolar decomposition of this ZFL re- 
sult in spin- weighted spherical harmonics. 

By computing the radiated momentum using Eq. (|4Tj) 
we find that it vanishes along the y— and z-axes, and 
that along the x-axis it is given by 

dP x pl 2 E 2 [v (15 - 13u 2 ) - 3 (v 4 - 6v 2 + 5) arctanhu] 
duj 3ttv 2 

(47) 

This result is in good agreement with the linear momen- 



tum radiated by point particles falling into Schwarzschild 
BHs, as we will see in section HVDI 



3. Zero-frequency limit 

Let us now consider the ZFL for generic values of the 
impact parameter. As bu — > we find once again that 
the energy spectrum is independent of the impact pa- 
rameter (as it was for the equal-mass collisions of sec- 
tion |nTAj. The leading-order expression of the energy in 
powers of fx/M is given by Eq. (|46p . reproducing Smarr's 
result for head-on collisions. Including higher powers of 
buj we get 



d 2 E 
dcudfl 



E 2 V ?v 2 



4-7T 2 (v sin cos </> — l) 2 
1 



v 2 (sin 2 cos 2 0-1) 



192 



v sin a cos 



( (372 - 8v 2 ) cos 26» + (6v 2 + 9) cos 40 + 2v 2 + 387) + 8 (2v 2 - 21) cos 20 



(20v 2 - 6) cos 40 - 6 (6v 2 + 35) 



{bu) 



O [(bu) 3 



r 



(48) 



{-4v (2v 2 + 3) sin 3 0(cos 20 + 3) cos 30 + 8 sin 2 cos 20 ( (3 - lOz; 2 ) cos 20 - 6v 2 + 9) + 



As in the equal-mass case, here too the radiation is sup- 
pressed along the x-axis, where the leading contribution 
is of order (btu) 4 : 



d 2 E 



dujdQ, 



E 2 y 2 {3-2v{v + Z)y 
576tt 2 



(M 



<D[{buf] . 
(49) 



C. Generality of the model 



The most important result of our ZFL calculation for 
collisions with generic impact parameter is perhaps that 
the ZFL itself is independent of the impact parameter. 
One of the limitations of the present calculation is that 
the modeling of the collision is rather ad-hoc, especially 
in the specification of the nature of the constraining 
forces. It is natural to ask how the results would change 
if the constraining forces were modeled differently For 
example, in our toy model the final system consists of 
two particles bound in a circular orbit, so the radiation 
spectrum shows peaks typical of the radiation produced 
by rotating bodies. There is a chance that the divergence 
at harmonics of the rotational frequency of the final sys- 
tem could contaminate the low-frequency behavior of the 
spectrum. 

To investigate this possibility, instead of considering 
the collision of two point particles, we studied a point 



particle colliding with a special extended matter distribu- 
tion: specifically, we considered an infinitely thin, slowly 
rotating, uniform disk. For brevity we do not report de- 
tails of this calculation here. Our main finding is that, 
if the disk is initially slowly rotating (so that after the 
collision the system is at rest), the ZFL is the same as in 
the case of two colliding particles. This is by no means 
a proof that the ZFL is completely independent of the 
way one models the system. It is however a hint that 
(as physical intuition would suggest) the ZFL should de- 
pend only on the asymptotic momenta of the colliding 
particles. 



IV. ULTRARELATIVISTIC INFALL OF POINT 
PARTICLES 

The formalism discussed in section IIIII is a flat-space 
approximation valid for the low-frequency part of the en- 
ergy spectrum and of the gravitational waveforms. In this 
section we compute the radiation from the linearized field 
equations in the curved background of a Schwarzschild 
BH. This is an accurate description at all frequencies in 
the limit where one of the binary components is much 
more massive than the other. 

There is extensive literature on the gravitational radia- 
tion emitted by p articles following geodesies in BH back- 
grounds (see |45j and Appendix C of [46| for summaries). 
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Here we present an incomplete overview of this literature. 
Davis et al. first studied the radiation emitted by par- 
ticles falling radially from rest into a Schwarzschild BH 
[47j and the synchrotron radiation emitted by particles 
in circular orbits lip (see also [4!| and Detweiler's con- 
tribution to Ref. [50(| ) . Oohara and Nakamura computed 
the energy, angular momentum and linear momentum ra- 
diated by particles falling from rest with generic angular 
momentum into Schwarzschild BHs (Hll HH . This work 
was later generalized to particles on scattering orbits in 
Schwarzschild, starting either from rest [53j or with finite 
energy at infinity [54j . Radial infalls into a nonrotating 
BH with finite energy were considered in Refs. [22l - |2~4j ] . 
and the ultrarelativistic limit was compared with Smarr's 
ZFL in Refs. [25T - [27j (in four dimensions) and [HI (in di- 
mensions D > 4). Oohara, Kojima and Nakamura stud- 
ied orbits plunging [jfj], [55l - l57j and scattering [58[ from 
rest in the case of rotating (Kerr) BHs. More recent stud- 
ies focused on the threshold of immediate merger using 
the geodesic analogy [Ta - tltij . 

Perhaps because of the limited astrophysical relevance 
of infalls with finite energy at infinity, to our knowledge 
there is no detailed study of the radiation emitted by 
point particles falling with generic energy and impact 
parameter into Schwarzschild BHs. One purpose of this 
section is to fill this surprising gap in the literature. A 
complementary study of scattering orbits with generic 
energy can be found in Ref. [54j . The generalization of 
this study to particles falling with arbitrary energy and 
impact parameter into Kerr BHs is in preparation. 

The radiation can be determined from the knowledge 
of the Sasaki-Nakamura wave function Xi m , which (in 
the frequency domain) can be written in the form 



L l in 



dr 2 



A 



u 2 - — (1(1 + l)r - 6M) 



Xlm — Sim ■ (50) 



Here (Z, m) are (tensor) spherical harmonic indices result- 
ing from a separation of the angular variables, u> is the 
Fourier frequency of the perturbation and A = r(r—2M). 
The boundary conditions dictate that we should have 
outgoing waves at infinity and ingoing waves at the BH 
horizon: 



Xi m — 



^lm e 
vout Auir* 
A lm e 



7'* —> 



-00 . 
00 . 



(51) 



The source term Si m in the Sasaki-Nakamura equa- 
tion (|50|) is determined by the point-particle trajectory. 
Without loss of generality we assume the trajectory to 
be an "equatorial" (6 = n/ '2) timelike geodesic in the 
Schwarzschild background, parametrized by 



2 62 



R 2 R 
R 2 



R 2 (E 2 



1) 



f(R)L 2 z 



f(R)f = E, 



2MR, (52) 
(53) 



where f(R) = 1 — 2M/R and dots stand for derivatives 
with respect to proper time t. As usual we denote by L z 
the orbital angular momentum of the particle along the 




FIG. 9: Trajectories for different values of L z /L CI it (as indi- 
cated in the legend) and E = 3. The black circle of radius 2 
marks the location of the horizon. 



z-axis and by E the particle's energy at infinity per unit 
mass /x, so E = 1 corresponds to an infall from rest. We 
remind the reader that the impact parameter b is related 
to L z by b = L z /y/ E 2 — 1. The geodesic equations can be 
integrated numerically for chosen values of E and L z . For 
L z > L cr j t = L CI it(E) (see section Hi A|) the particle does 
not plunge, but rather scatters to infinity. We restrict our 
discussion to the case L < L cl - it . Examples of plunging 
orbits for E = 3 and different angular momenta (i.e., 
different impact parameters) are shown in Fig. [3] 

Given a particle trajectory, the numerical calculation 
of the gravitational radiation emitted by the system in- 
volves the solution of Eq. (|50[) , which can be obtained by 
the standard Green's function technique (more details on 
the formalism can be found in Appendix [A]). Omitting 
for simplicity the indices (I , m) , we first define two inde- 
pendent solutions , -X^ut with boundary conditions 



*i 0) = 



r* — > — oo , 



Aw 



3 — iur* 



X, 



(0) 



, B in e~™ r * 
r* -» +oo , 



+oo 

— oo 



(54a) 
"(54b) 



Then the solution of the inhomogencous Sasaki- 
Nakamura equation (f50|) is given by 



Xu 



1 

W 



X 



(0) 



SimX^dr 



(55) 

where W = 2iujA ln is the Wronskian. Asymptotically for 
r* — > oo the amplitude of the wave function is 



1 

TP 



S lm X^dr 



(56) 



Then the radiated energy, angular momentum and linear 
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momentum are given by the multipolar sums 

dEi m 



did 

m dE lm 



Irn 

poo 

r " - L *■ s 

/"OG 

Px +iPy = / ^ 16w 2 V 



Hm^-l , m ^l m+ i T 



(57) 
(58) 
(59) 



Pirn \ A l ,m A l+l ,m+l ■'N .-m^M+l ,-m-l J 



with 

(IE; 



Pin 



16w 2 |A 



out 1 2 



(60) 



Qh 



(I + 3)(l - + to + 2)(l + m + l) 
(2Z + 3)(2Z + + l) 2 

y/(l-m)(l + m + l) 
1(1 + 1) 



1/2 



Note that the radiated momentum in the ^-direction 
P z = by symmetry. Details on the derivation of the 
source term and asymptotic expansions of the wave func- 
tions which are necessary to improve the numerical ac- 
curacy of the Wronskian are given in Appendix [X] 

We integrate all differential equations in CH — h using 
the adaptive stepsize integrator StepperDopr5 [59| . 
Schematically, the integration consists of the following 
steps: (i) integrate the first independent solution of the 
homogeneous SN equation ([5U|) with the boundary condi- 
tions (|54a[) from = 2M (1 + Sr) outwards (typically we 
choose 5r = 10~ 4 ); (ii) integrate the second independent 
solution of the homogeneous equation with boundary 



conditions (|54b[) from = /ui inwards, where typi- 
cally we choose 10 4 ; (hi) compute the Wron- 
skian at the large but finite radius rSS' , using Eq. (|A29[) 
for increased accuracy; (iv) integrate the geodesies with 
given orbital parameters and at the same time compute 
the source term using Eqs. (|A20[) and (|A23|) ; (v) out- 
put the "in" solution xj£' and the source term Si m on 
a grid of n = 1.6 x 10 5 colocation points, and use a 
Gauss-Legendre spectral integrator [59| to compute the 
convolution of the homogeneous solutions with the source 
term entering the expression for the outgoing amplitude 
X°^, Eq. (|5d| ; (vi) sum over multipoles to get the total 
radiated energy (|57[) , the angular momentum (|58[) and 
the linear momentum in the x— and y-directions (1591) . 

In the remainder of this section we summarize the re- 
sults obtained by this procedure. 



A. Point-particle spectra 

We performed an extensive set of simulations, select- 
ing seven values of the normalized particle energies (E = 



1, 1.5, 3, 5, 10, 20, 100). For each value of E we consid- 
ered nine different particle angular momenta (L z /L cr i t = 
0, 0.2, 0.5, 0.75, 0.9, 0.95, 0.99, 0.999, 0.9999), for a to- 
tal of 63 different configurations. We also ran a few more 
cases to validate our code against the results of Ref. [5l| . 
Our results are in good visual agreement with their plots. 
For each of the 63 simulations we computed all multipo- 
lar components of the radiation up to I — l max = 6 for 
10~ 2 < lj < 1.5 in steps of 8uj = 10~ 2 ; in a few selected 
cases we ran the code up to I = Z max = 8 to check the 
convergence of the results against truncation of the mul- 
tipolar sum. We also verified that for head-on collisions 
with E > 1 our results are in agreement (within three 
decimal places) with Table II of Ref. [24} . 

For each value of I, most of the radiation is typically 
emitted in the / = m component. In Fig. [TU] we show 
the I = m component of the energy spectra for particles 
falling from rest (E = 1) and for kinetic-energy domi- 
nated infalls (E = 3). In each of these two cases we plot 
the spectra for three selected values of the angular mo- 
mentum: L z /L CI it — (a head-on infall), L z jL a n =0.9 
and L z /L CI it = 0.9999 (a near-critical infall). The par- 
ticle trajectories corresponding to these values of L z are 
shown in Fig. |H] for E = 3. 

Not surprisingly, ultrarclativistic infalls radiate much 
more energy for a given value of L z /L cr i t . It is also ap- 
parent from Fig.[l0]that the energy output increases (and 
higher multipoles become relatively more important) as 
L z /L CT i t grows at fixed particle energy E. For kinetic- 
energy dominated infalls the energy spectrum approaches 
a nonzero constant as Mlu — > 0. The spectrum is flat 
for small frequencies in the head-on limit (25l - [28j . but 
the slope of the spectrum for Mco <C 1 is nonzero when 
the infall is nonradial. In all cases the energy spectrum 
decays exponentially at frequencies tu > Wq NM , where 
w qnm i s the fundamental Schwarzschild QNM frequency 
for the given multipole [13, HH . 

It is instructive to look at the multipolar structure of 
the radiation as we vary m for fixed I. Figure im shows 
spectra for the dominant quadrupolar component (I = 2) 
and all allowed values of m (|m| < I) for nonrelativis- 
tic infalls (E = 1) and kinetic energy-dominates infalls 
(E = 3) with L z /L crit = 0, 0.9, 0.9999. In the head-on 
case the (I + m)-odd components vanish, and compo- 
nents with the same |m| are exactly equal to each other 
because of symmetry. When L z the particle mo- 
tion breaks this degeneracy: the odd-m components are 
no longer zero, but they are still suppressed in the ZFL. 
The I = m (I = —to) components of the spectrum emerge 
from the ZFL with a positive (negative) slope, respec- 
tively. For large energies this slope is independent of the 
particle energy and proportional to the impact parame- 
ter, in qualitative agreement with the predictions of the 
toy model of section fTTTl 

For intermediate frequencies, the ZFL model of sec- 
tion Mil would predict a resonance in the spectrum at 
the rotational frequency of the rod, and a l/(u> — ujo) de- 
pendence near resonance. The perturbative spectra (as 
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FIG. 10: Spectra for I = m = 2, . . . , 6 for nonrelativistic infalls (E = 1, left column) and kinetic energy-dominates infalls (E — 3, 
right column). The top row refers to a radial infall, the middle (bottom) row to an infall with L/L cr i t = 0.9 (L/L C rit = 0.9999). 



well as the NR spectra shown in section [V] below) are 
instead characterized by an exponential decay for fre- 
quencies larger than the fundamental QNM frequency of 
the given multipolc. This difference can be attributed 
to the fundamentally different nature of the final state. 
In the toy model the final state consists of two parti- 
cles attached to a massless rod, with rotation frequency 
determined by angular momentum conservation (and an 
ad-hoc functional form of the stresses). In NR simula- 
tion and point-particle infalls, the post-plunge dynamics 



is dominated by the QNMs of the final BH. 

The m-even components all tend to the ZFL as Mui — > 
0, while the m-odd components are suppressed in the 
same limit. The dominant (/ = m) components of the 
spectra have a maximum corresponding to the QNM fre- 
quency and decay (roughly) exponentially for cj > Wq NM - 
In agreement with the toy model discussed in section Hill 
our pcrturbative results indicate that the ZFL of the 
spectra depends very weakly on the impact parameter, 
which (we recall again) is related to L z by Eq. @. This 
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FIG. 11: Spectra for I = 2 and |m| < I for nonrelativistic infalls (E — 1, left column) and kinetic energy-dominates infalls 
(E — 3, right column). The top row refers to a radial infall, the bottom row to an infall with L/L cr j t = 0.9. 



is shown very clearly in Table |TJ 

A close inspection of the right panels of Fig. [TT] shows 
that the ZFL of the spectrum for to = is not the same as 
for m = ±2. In fact, since the ZFL is almost independent 
of L z , one can use the head-on results to predict the rela- 
tive ratio between modes with m = and |m| = 2. In the 
head-on case, in a frame where the z-axis is aligned with 
the collision axis, only to = modes would contribute to 
the radiation. Once we rotate this "natural" coordinate 
system to the reference frame used in this paper, where 



a radial infall occurs along the x-axis (see Fig. [5]), for 
small L z the spectrum for (say) m = 2 will be related to 
the spectrum for to = by [60(: 

lim dE22/duJ = 2 . (61) 
u^o dE 20 /doj 2 v ; 

Our results for both equal-mass (see Fig. IT51 below) and 
extreme-mass ratio head-on collisions (Fig. Ill[) are in 
very good agreement with this prediction. 
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TABLE I: The numerically computed ZFL of the spectrum for 
different values of E is compared to the analytical prediction 
from a multipolar decomposition of the ZFL formula (|39|) . 
explained in Appendix [B] The agreement is remarkable. 



L z /L crit E = 1.5 E = 3 E = 10 E = 100 

0.000 0.0160 0.0481 0.0644 0.0663 

0.500 0.0160 0.0480 0.0643 0.0662 

0.750 0.0159 0.0480 0.0642 0.0662 

0.900 0.0159 0.0480 0.0642 0.0661 

0.950 0.0159 0.0479 0.0641 0.0661 

0.990 0.0159 0.0479 0.0641 0.0661 

0.999 0.0159 0.0479 0.0641 0.0661 

ZFL: 0.0158 0.0481 0.0644 0.0663 




FIG. 12: Infall from rest: spectra for I — m = 2 as L z /L CT i t — > 
1. In the inset: logarithmic divergence of the total radiated 
energy in the same limit. 

The energy spectra with L z /i cr ; t = 0.9, 0.9999 display 
secondary peaks which are related to the orbital motion 
of the particle. Radiation of orbital nature becomes more 
and more important as L z grows. Indeed, the nature of 
the spectra changes quite significantly as L z /L ct n — > 1. 
In Fig. [12] we show the I = m = 2 component of the 
spectrum for infalls from rest as we fine-tune the angular 
momentum to the critical value for a plunge: a very dis- 
tinctive "bump" appears at a frequency which is slightly 
lower than the QNM frequency, significantly enhancing 
the radiated energy. The location of this "bump" corre- 
sponds to (twice) the orbital frequency of the particle at 
the marginally bound orbit, i.e. w = 2f2 m b = (4M) _1 . 
For infalls from rest, a particle with L = L CT i t will orbit 
around the marginally bound circular geodesic. The rea- 
son for this local maximum in the spectrum is that when 
L z — >■ i cr ;t the particle orbits a large number of times 
close to the circular marginally bound orbit at r = AM. 

Gravitational radiation significantly affects near- 
critical geodesies for purely kinematical reasons. As 
L z /L CI it — > 1 the particle circles an infinite number of 



times around the marginally bound circular orbit with 
radius r = r m b, taking an infinite amount of proper time 
to reach the horizon. The proximity of the orbit to criti- 
cality is conveniently described by a small dimensionless 
parameter 

6=1- — . (62) 

-tJcrit 

As S — >• the particle hovers at a circular geodesic N ~ 
— _fficrit*crit log (kS)) I (T\/2V^f F ) times before plunging 
[61j . In our case V e s denotes the Schwarzschild effective 
radial potential, so that r 5 V^ = -24ML 2 +6rL 2 -4Mr 2 ; 
the angular velocity fi = dcj)/dt, dots stand for derivatives 
with respect to proper time and A; is a constant. All 
quantities are evaluated at the critical circular geodesic 
with L = L cr it, and radius r = r cr it- When E = 1 this 
circular geodesic corresponds to the marginally bound 
orbit, located at r m b = AM, MQ m b = 8 _1 and 

N ^=log<5. (63) 

In the ultrarelativistic limit E — > oo the critical circu- 
lar geodesic is located at the light ring r crit = 3M, the 
corresponding orbital frequency Mfi crit = (3-\/3) -1 and 

/V ~--L log (25) . (64) 

The orbital frequency at the light ring is intimately re- 
lated with the cikonal (long-wavelength) approximation 
of the fundamental QNM frequency of a BH [4l|, l62h - 
If54l |. This implies that ultrarelativistic infalls with a 
near-critical impact parameter are in a sense the most 
"natural" and efficient process to resonantly excite the 
dynamics of a BH. The proper oscillation modes of a 
Schwarzschild BH cannot be excited by particles on sta- 
ble circular orbits (i.e. particles with orbital radii r > 6M 
in Schwarzschild coordinates), but near-critical ultrarela- 
tivistic infalls are such that the orbital "bump" visible in 
Figure [T2] moves just slightly to the right to overlap with 
the "knee" due to quasinormal ringing. So ultrarelativis- 
tic infalls have just the right orbital frequency to excite 
BH oscillations. Resonant gravitational wave scattering 
explains the huge increase in radiated energy that can be 
observed in the bottom right panels of Figs. [10] and [TT] 
(see also (65l-[67|). Since QNMs are essentially perturba- 
tions of circular orbits at the light ring, we can expect 
that these conclusions will still apply to Kerr and even 
Kerr-Newman BHs [63|,[64[. 

For near-critical orbits, the bottom panels of Fig. ITTI 
show that the peak in the (I = 2, m = 2) component is 
located at twice the frequency corresponding to the peak 
in the (1 = 2, m = 1) component. The source term for a 
particle in circular orbit with frequency f2 typically con- 
tains a term proportional to 8(uj — mil); this lends further 
support to the "orbital" nature of the radiation enhance- 
ment. In fact, the inset of Fig. [T^] shows that the total 
radiated energy in the limit L z — > L cr it scales logarithmi- 
cally with 5, and hence linearly with the number of orbits 
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N, as expected for orbital radiation in the nonresonant 
case. In fact, in the next subsection we will show that our 
numerics are in quantitative agreement with the energy 
output of a particle in circular orbit at the marginally 
bound geodesic. 



B. Energy distribution 

We now turn our attention to the total integrated en- 
ergy. To simplify the discussion we start by revisiting 
and extending the analysis of infalls from rest (E = 1) 
first carried out by Oohara and Nakamura [5JJ . This 
reanalysis is useful both as a code check and to stress 
some important characteristics of the radiation in the 
near-critical limit L z /L cv n — > 1. Then we generalize our 
findings to infalls with arbitrary energy. 



1. Inf all from rest (E — 1) 

Any numerical calculation is necessarily limited to a 
finite number of multipolcs. In this paper we computed 
multipolar components of the radiation up to I = / max = 
6, but to get the total radiated energy, angular momen- 
tum and linear momentum we should in principle com- 
pute the infinite sums of Eqs. ([57jl. (|58]l and (HI]). This 
requires fitting the numerical results for I < l max by some 
analytical formula, and extrapolating this formula to es- 
timate the contribution from multipolcs with I > Z max . 
For infalls from rest, our results are consistent with the 
functional dependence proposed by Oohara and Naka- 
mura 1511: 



M F 



rad 



a E e 



-bm i 



(65) 



TABLE II: Fitting coefficients in Eqs. (JHSJ and (J72J). For 
each value of L z /L CI i t the first line refers to a fit including 
all multipoles, the second line to a fit dropping the I = 2 
multipole. A number such as 1.042(-2) means 1.042 x 1CT 2 . 



L Z J LcTlt 




b E 




aj 


bj 


l/^J rad 


0.000 





.45 


1.99 


1.04(-2) 


0.00 


0.00 


0.00 







.38 


1.95 


1.04(-2) 


0.00 


0.00 


0.00 


0.200 





.30 


1.66 


1.46(-2) 


1.71 


1.66 


8.20(-2) 







.23 


1.61 


1.46(-2) 


1.29 


1.61 


8.21(-2) 


0.500 





.31 


1.18 


4.50(-2) 


2.51 


1.24 


3.22(-l) 







.24 


1.13 


4.50(-2) 


1.75 


1.16 


3.23(-l) 


0.750 





.47 


0.92 


1.32(-1) 


3.60 


0.96 


8.98(-l) 







.37 


0.87 


2.32(-l) 


2.58 


0.90 


8.99(-l) 


0.900 





.76 


0.79 


2.92(-l) 


5.59 


0.83 


1.94 







.61 


0.75 


2.93(-l) 


4.14 


0.77 


1.95 


0.950 


1 


.02 


0.76 


4.27(-l) 


7.53 


0.80 


2.85 







.82 


0.72 


4.28(-l) 


5.63 


0.75 


2.86 


0.990 


1 


.76 


0.77 


7.30(-l) 


13.09 


0.80 


4.99 




1 


.36 


0.72 


7.32(-l) 


9.68 


0.74 


5.01 



A fit of the data yields the coefficients listed in Ta- 
ble HU The fitting coefficients listed in the first row are 
obtained by fitting all data (2 < I < 6). The lowest mul- 
tipole I = 2 is usually an outlier in this fit. As a rough 
check of the accuracy of the extrapolation, we repeat 
the fit considering only numerical data with 3 < Z < 6; 
this yields the coefficients listed in the second row. The 
difference between the total energies obtained by these 
two procedures can be seen as a very rough estimate of 
the error involved in the extrapolation. It is quite clear 
from the table that this error increases as L z /L cr ; t — > 1. 
This is of course a lower limit on the overall error in the 
computed energy, because it does not take into account 
numerical errors in the data, systematic errors coming 
from the (somewhat arbitrary) choice of the fitting func- 
tion, and inaccuracies in the fit itself. In any event, our 
fits are in good (but not perfect) agreement with the en- 
tries in Table I of Oohara and Nakamura [5l|. Oohara 
and Nakamura do not specify the range in I used for their 
fits, so it is hard to say if the slight difference in fitting 
coefficients is due to small differences in the numerics or 
to the fitting procedure itself. 

A fit to our perturbative results close to the critical 
angular momentum for capture yields 



prad = _o.0375 - 0.172 log 5. 



(66) 



where S was defined in Eq. (|62j). 

The analysis of the previous section reveals that close 
to the critical value of the angular momentum (i.e. as 
5 — > 0) the total radiated energy E lad should scale as 



— E™ d 



k + E N=1 N, 



(67) 



where Em=i is the energy radiated in one revolution 
close to the circular marginally bound geodesic, and k 
is an undetermined constant. Early calculations by De- 
tweiler [H, [Ho| (that we confirmed using a FORTRAN 
code discussed, for example, in Ref. [631) show that 
En = i ~ 0.65 for circular orbits with r cr i t = 4M in 
Schwarzschild. Using Eq. ([6"3"T) one gets the independent 
estimate 



M 

7- 



;E 



rad 



k- 0.15 log 5, 



(68) 



which is in good agreement with the fit of Eq. (166 



2. Ultrarelativistic mfall 

For ultrarelativistic collisions the scaling with I is not 
exponential, as illustrated in Fig. [13] In fact, we find 
that a power law of the form 



ad 



ME] 



c E l 



-cLe 



(69) 
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FIG. 13: Multipolar energy distribution for Ei/E 2 for E = 1 (left) and E = 100 (right). For large E the radiated energy scales 
like E 2 , not like jj 2 . Higher multipoles contribute (relatively) more when E is large. The scaling with I is well approximated 
by an exponential for collisions from rest and by a power law for relativistic collisions. 
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FIG. 14: Total energy radiated M/(fiE) 2 E lad rescaled by the 
particle energy squared, as a function of L z /L CI i t . In extrap- 
olating for / > 6, the fitting coefficients in Eqs. (|65|) and (|69[) 
were obtained either including (empty symbols) or excluding 
(filled symbols) the I — 2 data. The difference is only visible 
in the limit L z /L C rit — > 1. 



describes well our numerical results, where ce and d,E are 
constants. A least-squares fit yields the values of ce, cLe 
and the total radiated energy E lad listed in Table Mil 

For L z = and large E we get E" &d J(nE) 2 = 0.26, 
in agreement with the results of Refs. [25|, [28|. Close to 
the critical angular momentum (L z — > L C rit) the fit yields 
Ei ~ I /I. As explained in section |IVB[ in the limit E — s- 
oo the radiation is dominated by energy emitted at the 
light ring (the marginally bound circular geodesic located 
at r = 3M). Our fit is perfectly consistent with the 
classic study of synchrotron radiation by Davis et al. [48| , 
who found precisely a l/l dependence for the multipolar 
dependence of the radiation emitted by a particle orbiting 
at the circular null geodesic. 



For L z /L CT it < 0.95 we get the following fits for the 
total radiated energy as a function of the angular mo- 
mentum (see Fig. [T4|) : 



ME 



rad 



(HE)* 
ME™ d 

ME lad 



0.0685 exp 
0.145 exp 
0.294 exp 



3.241 



2.778 



2.176 



^crit 

L z 

L C rit 



,(-5 = 1.5), 
(£7 = 3), (70) 
(E = 100) . 



These results should be compared with the fit of the en- 
ergy computed from NR simulations of equal-mass BHs 



E r 



ad 



M 2 

ADM 



0.0195 exp 



2.632 



L z 



(E = 1.5) 



(71) 

In the fit to NR results (which have a numerical error of 
about 5%) the radiated energy is normalized to Madm 
and we estimate L z /L cr i t to be given by 6/6 C rit @- Un- 
fortunately, close to L cr it it is extremely difficult to get 
accurate estimates for £7 ad , because higher multipoles 
make an important contribution to the total radiation. 



C. Angular momentum 

1. In] all from rest (E = 1) 

Let us consider the angular momentum carried by the 
radiation. For slow motion, the multipolar decomposi- 
tion of the radiated angular momentum is consistent with 
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TABLE III: Fitting coefficients in Eq. (|69[). For each value of L z /L CT n the first line refers to a fit keeping all multipoles, the 
second line to a fit dropping the I — 2 multipole. 
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FIG. 15: Ratio M£ rad / Jrad for different values of E. The 
meaning of filled and empty symbols is the same as in Fig. 1141 
Oohara and Nakamura [5l|] find that this ratio is ~ 0.15±0.01 
for 1 < L z < 3.9 when E — 1 and L cr it = 4, but we expect 
that £ rad /J rad ~ (Mfi 8cat ) as L z /L clit -> 1 (see text). 



the exponential dependence 



I 1 



rad 



aje 



-bj I 



(72) 



where the coefficients aj, bj are given in Table HU 

As shown in Fig. [15l our results for infalls from rest are 
consistent with Oohara and Nakamura's [5l| suggestion 
that for 0.25 < L z /£ C rit S 0-99 the radiated angular 
momentum should obey the approximate relation 



ME 1 



rad 



J 



rad 



(0.15 ±0.01) 



(73) 



However, this approximate proportionality relation must 
break down as L z — > L CIlt , because for aparticle in circu- 
lar motion E Tad /J Iad ~ (Mfl scat ) [H,^. For example, 



FIG. 16: Total radiated angular momentum (obtained by ex- 
trapolation) as a function of L z / L CI i%. The meaning of filled 
and empty symbols is the same as in Fig. 1141 



for E = 1 we should find that E Tad /J Tad -> 1/8 = 0.125, 
and as E -> oo we expect that £ rad /J rad -> l/(3\/3) ^ 
0.19245. Quite remarkably, and despite the sizeable un- 
certainties in extrapolating our numerical results to get 
the total radiated energy and angular momentum, Fig. 1151 
is consistent with these predictions. 



2. Ultrarelativistic mfall 

The total radiated angular momentum for generic infall 
energies and its multipolar decomposition are shown in 
Figs. lTBlfFTl For collisions with generic impact parameter, 
we find that the angular momentum emitted in a given 
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FIG. 17: Multipolar components of the radiated angular momentum Ji/E 2 for E = 1 (left) and E = 100 (right). For large E 
the radiated angular momentum scales like (fiE) 2 . Higher multipoles contribute (relatively) more when E is large. The scaling 
with I is logarithmic for collisions from rest, but it acquires corrections for rclativistic collisions. By symmetry, no angular 
momentum is radiated when L z — 0. 



multipole I has a power-law dependence of the form 

Trad 

Ttw = cjl ~"' (74) 

where cj and dj are constants. The values of cj, dj and 
of the total radiated angular momentum J rad obtained 
by extrapolation are listed in Table HVl The angular mo- 
mentum emitted in each multipole I < Z max is shown in 
Fig. [T7] Using Eq. ([74)) we can extrapolate to get the 
total radiated angular momentum. Fig. [TS] shows that 
the ratio J rad /(M£' rad ) for kinetic -energy dominated in- 
falls is only (roughly) constant in the intermediate regime 
0.5 < L z /L cr it < 0.95. The following expressions provide 
good fits to the extrapolated values: 



:atl 



jrad 
jrad 



= 0.0142 



0.0676 



1 - exp 3.170 



1 - exp 2.585 



£crit 

L, 



(E =1.5). 



(£7 = 3) 



0.42 



1 - exp 1.836 



L, 



, (£7 = 100) 



(75) 



These fits can be compared with the corresponding fit 



from NR simulations of equal-mass BHs 
(within a numerical error of about 5%) 



which yields 



rad 



M 2 

ADM 



0.0024 



1 - exp 2.928 



n 2 



(E 



1.5). 
(76) 



D. Linear momentum 



multipoles contribute in building up the total linear mo- 
mentum, and also to "visually" test the convergence of 
the sum in Eq. (|5"9"|) . we plot in different linestyles the 
spectra obtained by summing over the lowest ^ max multi- 
poles with i max = 2 , . . . ,6. Strictly speaking, if we sum 
Eq. (|59j) up to Z max = 6 the last term of the sum will 
be inconsistent, because it involves amplitudes such as 
^z+i.m+i, which we have not computed and we set to 
zero for lack of a better alternative. Note that this "mul- 
tipolar coupling" does not occur in the calculation of the 
energy and angular momentum. The inconsistent trun- 
cation means that we should only consider these plots as 
representative of the real convergence properties of the 
sum for Z max < 5. Furthermore, the spectrum obtained 
when we truncate at Z max = 6 cannot be trusted in fitting 
numerical results to extrapolate the sum to infinity. For 
this reason (and also for the highly oscillatory buildup of 
linear momentum which is evident from Fig. I18j) the ex- 
trapolation of numerical results to get the total radiated 
linear momentum is quite sensitive to the relatively low 
number of multipoles that we are using. 



1. Inf all from rest (E = 1) 

Despite these caveats, in the case of low-energy col- 
lisions (and in particular for infalls from rest) the total 
linear momentum converges reasonably fast with I. Sum- 
ming over 2 < I < 6 is typically enough to warrant an 
accuracy of 10% or less in the total linear momentum. 

Nevertheless we have used extrapolation, fitting our 
numerical results by an exponential: 



/max = " _ ttPe 



bp l n 



(77) 



Representative spectra of the radiated linear momen- 
tum are shown in Fig. 1181 To illustrate how different 



The fitting coefficients and the total radiated momentum 
obtained by extrapolation are listed in Table [V] 
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TABLE IV: Fitting coefficients in Eq. (|74|) . For each value of L z /L CI it the first line refers to a fit keeping all multipoles, the 
second line to a fit dropping the I — 2 multipole. 



E = 1.5 E = 3 E = 10 E = 100 

LJL^t cj dj l/(fiE) 2 J™ d cj dj l/(fiE) 2 J™ d cj dj l/( M g) 2 J rad cj dj l/( M g)V rad 






000 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 





00 


0.00 


0.00 





00 


0.00 






0.00 


0.00 


0.00 


0.00 


0.00 


0.00 


0.00 





00 


0.00 


0.00 





00 


0.00 





200 


3.79 


4.60 


0.16 


1.30 


2.71 


0.35 


0.82 


1 


98 


0.54 


0.78 


1 


91 


0.58 






11.49 


5.32 


0.16 


1.89 


2.95 


0.34 


0.83 


1 


99 


0.54 


0.76 


1 


89 


0.58 





500 


5.22 


3.58 


0.58 


3.11 


2.45 


1.10 


2.20 


1 


89 


1.68 


2.09 


1 


83 


1.75 






9.62 


3.98 


0.57 


4.15 


2.64 


1.08 


2.24 


1 


90 


1.67 


2.05 


1 


82 


1.76 





750 


5.55 


2.70 


1.50 


4.48 


2.09 


2.52 


3.62 


1 


72 


3.65 


3.42 


1 


69 


3.69 






7.60 


2.90 


1.47 


5.40 


2.21 


2.45 


3.72 


1 


74 


3.61 


3.38 


1 


68 


3.71 





900 


6.59 


2.15 


3.41 


5.76 


1.79 


5.11 


5.11 


1 


58 


6.88 


4.84 


1 


56 


6.70 






8.09 


2.28 


3.32 


6.46 


1.87 


4.94 


5.25 


1 


60 


6.77 


4.82 


1 


56 


6.72 





950 


7.77 


1.93 


5.51 


6.65 


1.63 


8.09 


6.11 


1 


48 


10.40 


5.83 


1 


47 


10.08 






9.43 


2.06 


5.29 


7.30 


1.69 


7.75 


6.33 


1 


50 


10.13 


5.89 


1 


48 


9.98 





990 


11.85 


1.71 


12.15 


9.02 


1.35 


21.65 


8.45 


1 


26 


27.81 


8.14 


1 


26 


27.54 






15.02 


1.86 


11.22 


9.96 


1.42 


19.52 


8.98 


1 


30 


25.35 


8.53 


1 


29 


25.55 



TABLE V: Fitting coefficients in Eq. ([77]). The first line refers 
to a fit including all multipoles, the second line to a fit drop- 
ping I — 2. 



L z 1 LcTlt 


ap 


bp 


A///i 2 P rad 


0.000 


7.03(-3) 


2.00 


8.33(-4) 




7.98(-3) 


2.05 


8.33(-4) 


0.200 


1.09(-2) 


1.56 


2.41 (-3) 




9.99(-3) 


1.53 


2.42(-3) 


0.500 


3.05(-2) 


1.13 


1.16(-2) 




2.58(-2) 


1.06 


1.16(-2) 


0.750 


6.54(-2) 


0.84 


3.37(-2) 




5.55(-2) 


0.78 


3.39(-2) 


0.900 


9.24(-2) 


0.67 


5.75(-2) 




8.13(-2) 


0.61 


5.81(-2) 


0.950 


8.85(-2) 


0.58 


6.20(-2) 




7.83(-2) 


0.52 


6.30(-2) 


0.990 


6.77(-2) 


0.72 


4.96(-2) 




4.03(-2) 


0.48 


5.16(-2) 



2. Ultrarelativistic infall 

Spectra of the linear momentum radiated by ultrarel- 
ativistic infalls are shown in Fig. [THJ For head-on colli- 
sions, in the limit of large boosts the radiated momentum 
is well described by 

i~-(o.ii-|&)f#, — p., 

For ultrarelativistic infalls the ZFL of the spectrum 
is nonvanishing and it depends very weakly on the im- 
pact parameter (we observed a similar trend when study- 
ing the ZFL of the spectrum of the radiated energy). 
The ZFL is very well approximated by Smarr's formula, 
Eq. (jTT)) . This is shown quantitatively in Table IVI1 



TABLE VI: The ZFL of the radiated linear momentum for 
head-on collisions, compared with the analytical result (|47[) . 







i 






L z 1 Lcrit 


E = 1.5 


E = 3 


E = 5 E = 10 


E = 100 


0.00 


0.0182 


0.100 


0.151 0.187 


0.203 


ZFL: 


0.0180 


0.0993 


0.150 0.189 


0.212 



We tried to fit our numerical results for high energy 
and generic impact parameters by a power law of the 
form P/ max = P rad — cp?~fj. For the reasons explained 
above, the errors associated with this extrapolation pro- 
cedure arc significant, and wc decided not to present the 
extrapolated values of _p rad . 

Qualitatively, in the ultrarelativistic case the total ra- 
diated linear momentum increases only mildly with the 
impact parameter: we find a maximum increase of a fac- 
tor ~ 3 relatively to the head-on case. This observation 
can be used to give a rough estimate of the maximum 
recoil that could result from gravitational- wave emission 
in the ultrarelativistic collision of nonspinning BHs. If 
we take fj,E/M = 1/10 we get a recoil velocity of or- 
der v ~ 10~ 3 for head-on collisions, and a maximum 
of around v ~ 3 x 10 -3 for grazing collisions. Restor- 
ing physical units, this would correspond to recoils in 
the range 300 — 900km/s. It will be interesting to ver- 
ify these estimates by NR simulations of ultrarelativistic, 
comparable-mass BH binaries. 



V. COMPARING DIFFERENT APPROACHES 

As discussed in Ref. Q (see also Fig. [I}, NR results 
for the energy spectra from the high-energy, head-on 
collision of two BHs are in very good agreement with 
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FIG. 18: Spectra of the linear momentum radiated for infalls from rest (E — 1) and ultrarelativistic infalls (E = 100). Different 
lines refer to different truncation choices in the sum of Eq. (|59|) . In the head-on case (top row) the only nonzero component of 
the linear momentum is P x ; for L/L cr i t = 0.9 (middle and bottom rows) both the P x and P v components are nonzero. To test 
convergence we sum Eq. (|59[) up to different values of i max , as indicated in the legend. 



the ZFL predictions. Calculations of the radiation from 
point particles falling into nonrotating BHs are in quan- 
titative agreement with ZFL predictions in the extreme- 
mass ratio limit (see for instance Table H| , and in qual- 
itative agreement with NR calculations. In this section 
we present a more extensive comparison between these 
different approaches. 



A. Head-on collisions 

For head-on collisions, the ZFL results of section IIIII 
are in qualitative and quantitative agreement with point- 
particle calculations. In particular, a multipolar decom- 
position of ZFL calculations yields a flat energy spectrum 
for all multipolar components (see Appendix |B|) . Point- 
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particle calculations and NR simulations suggest that the 
cutoff frequency for each multipole can be chosen to be 
the lowest QNM frequency for the given multipole. A 
multipolar analysis is therefore important to introduce 
a more "natural" and appealing cutoff frequency in the 
ZFL spectrum, in contrast with Smarr's original, quali- 
tative suggestions [20j |. 

In the eikonal limit, the real part of the fundamental 
QNM frequency for the l-th multipole Wq NM is related to 

the frequency M£l c = (3\/3) _1 of unstable circular null 
geodesies [UlSD]: 



QNM 



(79) 



Combining this approximation (which is surprisingly 
good even for low values of I) with ZFL calculations of 
the energy spectrum at low frequencies, we can estimate 
the total radiated energy as follows: 



E 



rati 



1=2 



x w, 



QNM 



(80) 



At high energies the ZFL assumes a particularly simple 
form. Using the multipolar decomposition results from 
Appendix [Bl we get, for the equal-mass case: 



E 



rati 



E 

even 1 



4(EM) 2 (2/ + l)(Z - 2)! 



4 log 2 



EM 



(1 + 2)1 
0.1698EM. 



n) 



Here we used EMVL C = (3A/3)" 1 , since the final BH 
has (approximately) mass EM. An extrapolation of re- 
sults from recent NR simulations Q predicts E lad ~ 
(0.14 ± 0.03)MadMj hi remarkably good agreement with 
this naive estimate. 

As another example, let us consider the extreme mass- 
ratio case in the high-energy limit. Using again results 
from Appendix [B] we get 



E 



rad 



y, 4£V (2/ 



l)(Z-2)! 



13 



7T 



9V3tt M 



(I 

0.265 



2)! 

p?E 2 
M 



in c 



(82) 



The extrapolation of point-particle results to E —> oo 
(see iEUll and Table EH predicts £ rad - 0.262fi 2 E 2 /M 
[251 ] , again in remarkable agreement with our simple ap- 
proximation. 



B. Collisions with finite impact parameter 

Unfortunately, an extension of our head-on estimates 
to the generic case of collisions with finite impact param- 
eter is not straightfoward, because the energy spectrum is 
no longer flat. However several generic features are com- 
mon to the ZFL, point-particle and full NR calculations. 



Fig.ll9lillustrates our point. In the left panel we show the 
energy spectra (rescaled by the total ADM mass) for NR 
simulations of equal- mass BH collisions with E = 1.5 and 
varying impact parameter . In the right panel we show 
spectra for point particles falling into a Schwarzschild 
BH of mass M with energy E = 1.5 and different im- 
pact parameters. The angular momentum of the particle 
normalized by the critical angular momentum (which for 
E = 1.5 is L cr it ~ 6.35M) is indicated in the legend. 
Strictly speaking, these two plots can only be compared 
from a qualitative point of view. In the equal-mass case 
the final BH is spinning (in fact, for large impact parame- 
ters it is rapidly spinning [(| ) , whereas our point-particle 
calculation considers nonspinning holes. Nevertheless, 
the spectra show qualitative agreement. In particular, 
the ZFL is independent of the impact parameter for both 
equal-mass and extreme-mass ratio collisions. For point 
particles falling into Schwarzschild BHs the cutoff fre- 
quency does not depend on m, but early calculations of 
nonrelativistic equatorial infalls into Kerr BHs j55l - [57| 
hint that this degeneracy in the cutoff frequency should 
be lifted when one considers the rotating case. 

For small, but finite frequencies, we find that the spec- 
trum of positive-m (ncgative-m) modes has positive (neg- 
ative) slope as u) — > 0. The ZFL-inspired calculation 
for generic impact parameter of section IIIII taught us 
that while the ZFL itself is a robust feature, the finite- 
frequency behavior is not, and it strongly depends on 
the modeling of constraining forces. Such arbitrariness is 
absent in the point-particle calculation. 

A fit to NR results for E — 1.5 and I = m = 2 yields 



1 dE 
~M 2 ~~d^ 



1 dE 
lvP~d^ 



0.15- 



-Moj. 



(83) 



w=0 Jj cr it 

In the point-particle limit, a fit of the I = m = 2 spectra 
with E = 1.5 yields 

1 dE 1 dE 



(HE) 2 duj (^lE) 2 duj 



0.65-^ Mui. 
E cr it 



(84) 



VI. OUTLOOK 

In this paper we have used a combination of ZFL cal- 
culations and pcrturbative techniques to study some of 
the main features emerging from ongoing NR simulations 
of ultrarelativistic BH collisions. Here we wish to point 
out possible extensions of our analysis. 

We are currently working on a perturbative analysis 
of ultrarelativistic infalls with generic impact parameter 
into Kerr BHs (see [5514571 ] for studies of infalls from rest). 
Another obvious generalization would be to extend ZFL 
and point-particle calculations to higher-dimensional 
spacetimes. A preliminary investigation of head-on, ul- 
trarelativistic infalls into D-dimensional (Schwarzschild- 
Tangherlini) BHs can be found in Ref. [28[. It would be 
interesting to extend that study to plunges with generic 
energy and impact parameter, and eventually also to ro- 
tating (Myers-Perry) BHs. 
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FIG. 19: Left: I — m = 2 component of the energy spectrum of the gravitational radiation emitted in the collision of two 
equal-mass BHs having speed /3 = v/c ~ 0.75 in the center-of-mass frame 0,01- All quantities are normalized to the ADM 
mass of the system Madm- Right: I = m = 2 component of the energy spectrum of the gravitational radiation emitted by 
point particles falling into Schwarzschild BHs of mass M with energy E = 1.5. 



As mentioned in the introduction, NR simulations have 
provided evidence of the existence of two critical thresh- 
olds, depending on the impact parameter b of the colli- 
sion: for b < b* (where b* is the threshold of immediate 
merger) the BHs merge within the first encounter; for 
b* < b < 6 sca t (where 6 scat is the scattering threshold) the 
binary does not merge immediately, but sufficient energy 
is radiated to put the binary into a bound state that 
eventually results in a merger; finally, for b > fegnat the 
BHs scatter producing bremsstrahlung radiation [6( . The 
emergence of these two different thresholds &* and & sca t 
can be explained in terms of a gravitational radiation- 
induced splitting of the scattering threshold. This split- 
ting can be described in terms of the so-called Mclnikov 
function, which is well-known from the theory of dynam- 
ical systems [7(| • A detailed analysis of this problem will 
be the topic of a future publication. 

In this paper we have not discussed the radiation 
emitted by ultrarelativistic encounters leading to scat- 
tering, rather than BH formation. Studies of gravi- 
tational bremsstrahlung started in the sixties with the 
work by Peters [7l|, and they employed several approx- 
imation schemes. Turner and Will studied the gravita- 
tional radiation emitted as a result of the scattering in 
a low- velocity (post- Newtonian) approximation [72h 74 1. 
D'Eath and collaborators developed an elegant approxi- 
mation to estimate the radiation in the ultrarelativistic 
case (75l480j . Matzner and Nutku developed a "method 
of virtual quanta" analogous to the Wcizsacker- Williams 
method of virtual quanta in electromagnetism, which is 
restricted to high-velocity, extreme-mass ratio encoun- 
ters [ll| . The most complete study of bremsstrahlung is 
perhaps a series of papers by Kovacs, Thorne and Crow- 
ley, who computed the radiation from compact objects of 
arbitrary mass ratio flying past each other at arbit rary 
velocities in the limit of large impact parameter [82h85] . 



We refer the reader to section VI of Rcf. [85J for a com- 
prehensive review of methods of computing gravitational 
bremsstrahlung radiation. 

In the eighties, various groups worked out the gravita- 
tional radiation from particles scatterin g ar ound BHs at 
nonrelativistic and relativistic velocities [53l.l54j (see also 
[4o| for a review) . It would be interesting to repeat their 
perturbative calculations paying special attention to ul- 
trarelativistic scattering and to the near-critical behav- 
ior. Moreover, a comparison of bremsstrahlung radiation 
as computed in NR against all these different approxima- 
tions would be a very interesting topic for future work. 

Yet another semianalytical approach that we have not 
considered in this paper, but that could certainly prove 
useful as a diagnostic of NR codes in D dimensions, 
is the close-limit approximation (see e.g. [HI). A de- 
tailed quantitative understanding of gravitational radi- 
ation from Z?-dimensional BH collisions will ultimately 
rely on the development of NR in higher dimensional 
spacetimes. Vario us g roups are making rapid progress 
in this direction IS-QjJ. 
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Appendix A: The Sasaki-Nakamura formalism 



Here and below we generally omit the subscripts (I, m) 
to simplify the notation. We also define A = r 2 f(r) = 
r 2 — 2Mr, p = —1/r and introduce the differential oper- 
ators 

vn 

C s = d e + — - + scot 6>, (A7) 

SlUtf 

J± = dr±i%. (A8) 



The gravitational radiation generated by point parti- 
cles in BH backgrounds is best described in the Sasaki- 
Nakamura formalism (87l - [89| . Here we summarize the 
equations describing the general infall of particles of rest 
mass /i with arbitrary energy (per unit rest mass) E. We 
consider the metric of a nonrotating BH in Schwarzschild 
coordinates: 



ds 2 = -f(r)dt 2 + fir^dr 2 + r 2 (d9 2 + sin 2 

(Al) 

where f(r) = 1 — 2M/r and M is the BH mass. In the 
Teukolsky formalism the perturbation equations can be 
reduced to a second-order differential equation for the 
Newman-Penrose scalar ip^ with a source term Tt- We 
can expand ^4 as 

dw^Rlmir) -2Y lm (n)e- iut , 

-°° lm 

(A2) 

and similarly for Tt- We denote by s Yi m {9) the spin- 
s weighted spherical harmonic, which can be expressed 
in terms of the well-known (scalar) spherical harmonics 

For a point particle falling into a BH along a 
geodesic the coordinates can be parametrized in terms 
of (say) the radial location of the particle: (t, r, 0) = 
(T(R), R, Q(i£)). The source term in Teukolsky's equa- 
tion is directly related to the energy-momentum tensor 
T^ v of the test particle of mass p: 

T^ = ^5(t-T(R))5 2 (n-m))^r^r, (A3) 
r z \r\ dr dr 

where an overdot denotes d/dr. In particular, 
Tt = 4 1 dQdtp- 6 (B' 2 lm + B* 2lm ) e-^^YUOy 



(A4) 



The quantities T nn = T^n^n", T nm = T^n^fh^ , and 
Tmfh = T liv fh^fh v are contractions of T^ v with the Kin- 
nersley tetrad "legs" 



1 



(/(r), 1,0,0) , 



(A9) 



(0,0,-r 2 ,ir 2 sin6>) , (A10) 



and an overbar denotes complex conjugation. The ex- 
plicit expressions are 



T »« = £2T^(f(r)i + r) 2 S(t-T(R))5 2 (n - Q(R)) , 

Tnm = (f(r)i + f) S(t-T(R))S 2 (Q-Q(R)), 

Tram = -~5(t~T(R))S 2 (n-n(R)). 
2 \r\ 



It is convenient to move the delta functions out of the 
operators by performing an integration by parts. To do 
that, note that given two functions A and B 



dttA£ s B 



dflB a 



sm( 



(s - l)cot0j A, 
(All) 



where we assumed that the integrals exist and we defined 



c+ = a, - 



sin 9 



s cot 9 . 



(A12) 



with 



From the gen eral properties of spin-weighted spherical 
harmonics [9(| we get 



B' 2lm = --p 9 C^[p- 4 C Q {p- 3 T nn ) 



2V2 



p 9 A 2 C-i[p~ 2 J + (A- V 4 ^)] ,(A5) 



B* 2 ' lm = --p 9 A a J + [p- 4 J+ (p" 1 ^)] 

+ ^p 9 A 2 J + [p- 2 £_x (A-V- 4 T„m)] .(A6) 



C+C+ - 2 Y lm = y/X(X + 2) Y t 

-^2 -2^m = 



(A13) 



y/(l + 2)(l-l) -iY lm , (A14) 



where A = (I — 1)(Z + 2). Introducing the operator L + = 
■A — h ioJ, where the tortoise coordinate r* is defined by 

^ = ^, we can write T T = T x + T 2 + T 3> with 



2G 



/. . Jl —imip+iuj t 
dndt(f?C-i\p-*Co (p- 3 T„ n )]) e ^ 4 -^_ 2 y ;m (0) = -ft- (f(r)i + ryC+ C+- 2 Y lm {6), 

= -JLA 2 {V'f r VA(A + 2)e iuT - im ^ Q Y lm {6 ) , (A15) 

t 2 = 4 y ^^e i " t - m ^_ 2 y im (0)[£_ 1 ( P -V + (A-V^ 4 r„„-O) + J+(p- 2 £-i(A-V~ 4 ^))] 

T 3 = ~y dfie^*- < '^_ 2 y Im (fl)AVJ'+(p" 4 J r +(p" 1 r fiWft )) =^i+ (x £ + (^^ efa,r ~ <m *)) -a?im(0o) • (A17) 



Here we introduced the advanced coordinate V = i + r* 
and we denoted radial derivatives with a prime (so V = 
dV/dr). 

In the Sasaki-Nakamura (SN) formalism [87j-t89| one 
introduces a new wave function Xi m {r) related to the 
radial Tcukolsky function of Eq. (|A2[) via 



Ri 



AL+ ( -L+ (rX lm ) 



(A18) 



The function X[ m satisfies the differential equation 
where the source term Si m is related to the Teukolsky 
source Tt by the relation 



AL+ I — ] L 



r 5 S lr 
A 



-T T . 



It is convenient to define a quantity 

„5 



which satisfies 



W lr . 



d 2 W h 



— «?, p lbjr 



r 

-— T T e l 
A 2 



(A19) 



(A20) 



(A21) 



(A22) 



dr 2 

By using the equality 

L+ (X im (r)e^ T ) = d r , (X lm (r)e iuV ) 
we can write Wi m = W% + W% + W3, where 

^ = v/A(A + 2)| (T/') 2 re^-™* o y, m (0 o ) , 

i/iVA^.(v"r a 0e to ' v - 4m *) _i^ m (fl ). (A23) 

-2^771(^0) • 



dr 2 
d 2 W 3 
dr 2 



Given any function F(r) we have F(r)V 'e wV = 
l/(iw)[( J F 1 e^ y )' - F'e iuV }. If we use this identity to in- 
tegrate by parts, rearrange terms and set f = —7, these 
equations reduce to Eqs. (A5)-(A8) in Ref. 51 1. 



1. Calculation of the Wronskian 



From a numerical perspective, the integration of 
Eq. (|50[) is performed via a standard Green's function so- 
lution. For improved accuracy, in the numerical integra- 
tions we use asymptotic expansions of the wave functions 
near the horizon and near infinity. Close to r = r + = 2M 

we set y = r ^ ; make the ansatz 



X lm (r) 



OO 



k=l 



1 + } akV" 



(A24) 



and substitute into the differential equation (|50p . Noting 
that / = y — y 2 + 0(y 3 ) we determine the leading-order 
and next-to-leading-order coefficients to be 



01 



«2 



A - 1 
1 - AiMuj ' 

3+ (A - l + 8iMu))ai 
4(1 - 2iMw) 



(A25) 



In the limit r — > 00 we have 



k=l 



fc=i 



(A26) 

The lowest-order coefficients in the series are found to be 



*(A + 2) A(A + 2) + 12iMu 
a = — ^ , c 2 = ^-9 i ( A27 ) 



2u 
-ci , d 2 



8uj 2 

UiMu - A(A + 2) 
8U 2 



(A28) 



One can invert these relations to get 
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[(2M - r)(2c 2 + cir) + ir 2 (c 2 + r(ci + r))u] X lm (r) - r (r - 2U)[c 2 + r(ci + r)]X,' m (r) 
(2M - r) [cid 2 + 2d 2 r + (di - ci)r 2 - c 2 (di + 2r)] + 2r[c 2 + r(a + r)][d 2 + r(di + r)]u 

[i (2M - r)(2d 2 + dir) + r 2 (d 2 + r(di + r))tj X lm (r) -ir(r- 2M)[d 2 + r(di + r)}X' lm (r) _. 

(2M - r) [-cicfe - 2d 2 r + (ci - di)r 2 + c 2 (di + 2r)\ + 2ir[c 2 + r(a + r)][d 2 + r(di + r)]u & 



where a prime stands for a derivative with respect to r 
and in our numerics all quantites are evaluated at large 
but finite r (typically at r = r x = r£) /u, with r& = 4x 
10 4 ). The Wronskian appearing in the Green's function 
solution of the inhomogeneous equation can be evaluated 
from the relation W = 2iu)A l ^L . 



Appendix B: Multipolar decomposition of ZFL 
spectra 

Here we discuss the extraction of multipolar compo- 
nents from the ZFL calculation. For consistency with the 
conventions we used in NR simulations [3, @ and in the 
point-particle infalls discussed in this paper, we assume 
that the polarization states have an angular dependence 
that can be decomposed in spin-weighted spherical har- 
monics: h + — ih x ~ ^2 him- 2 Yi m . Then (using the com- 
pleteness of spin-weighted spherical harmonics) we can 
extract the multipolar content of the energy spectrum 
by equating 



d 2 E 
dwdQ 



E 

V lm 



dE, 



1 1 II 



du) 



(Bl) 



where dEi m /dui are (as yet undetermined) functions of 
u>. From the orthonormality of spin- weighted spherical 
harmonics wc find 



dEi„ 
dui 



dQ 



d 2 E 
diodQ 



--iYu. 



(B2) 



We now determine analytically the multipolar content of 
the radiation in the simple case of ultrarelativistic head- 
on collisions. 



1. Equal mass collisions 

If we assume that the collision occurs along the z-axis, 
the radiation will only contain m = modes. For equal- 



mass, ultrarelativistic head-on collisions Eq. (|39[) implies 
dEio 4E 2 M 2 (2l + l)(l -2)\ 



du 



(1 + 2)1 



= 0, 



/ even 
/ odd. (B3) 



Odd multipoles do not contribute, as required by equa- 
torial symmetry. Summing over even multipoles we get 
dEjduj = E 2 M 2 /n, in agreement with Smarr's Eq. (2.20) 
when v — > c [2(|. This expression assumes that the colli- 
sion occurs along the z— axis, but the multipolar compo- 
nents in a general (rotated) frame can be found following 
the procedure discussed in Ref. [oOj . 



Extreme-mass ratio collisions 



From Eq. (gBJ}, for fi = Mi <C M 2 = M we get 



dEip _ 4E 2 fi 2 (2Z + — 2)! 

did 7T 



(1 + 2)1 



(B4) 



Summing over multipoles we get dE/du = 4/ (3tt)E 2 u 2 , 
in agreement with Smarr's Eq. (2.18) when v — > c j20j . 

We recall once more that Eq. (|B4[) is valid in a frame 
where the collision occurs along the z-axis, so only m = 
modes are present. The transformation to a general (ro- 
tated) frame is explained in Ref. [60[. In the coordinate 
system used to compute the radiation from point par- 
ticles falling into a Schwarzschild BH (section ITV]) , and 
focusing on the dominant (I = 2) components, we get: 



dE 20 
du> 

dE 22 
du 



5E 2 n 2 
24tt 

5E 2 fi 2 
16tt 



0.0663146.EV , 



0.0994718.EV 



(B5) 
(B6) 



This is in excellent agreement with the point-particle re- 
sults listed in Table fl] 
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